Tarea 5: Bayesian Inference Part II

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Fecha límite de entrega: 13/12/2023

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Primera Parte: Preguntas Teóricas
  5. Segunda Parte: Elaboración de Código

Objetivo

a la uuuuultima tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la ultima parte del curso, los cuales se enfocan principalmente en aplicar inferencia bayesiana para generar regresiones lineales y estudiar métodos de obtención de la posterior mas poderosos, como es MCMC. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Documentación:

Pregunta 1: Model Evaluation and Information Criteria

Explique cómo cross-validation, criterios de información y regularización ayudan a mitigar o medir los problemas de underfitting y overfitting.

Cross-Validation

La cross-validation divide los datos de entrenamiento en \(k\) subconjuntos, entrenando el modelo \(k\) veces, siempre dejando uno de los subconjuntos fuera. Esto logra que el modelo no se ajuste a unos datos en particular, sino que aprenda de mejor manera de todos los datos.

Criterios de Información

Los criterios de información permiten, de cierta forma, medir cuanto overfitting o underfitting tienen distintos modelos entrenados con los mismos datos. Esto permite comparar distintos modelos para ver cual es el que se comporta mejor. Aún así, esto no resuelve el problema de overfitting, solo lo mide, ya que estas medidas mejoran en general mientras más complejo (y por lo tanto, más ajustado) es el modelo a los datos.

Regularización

La regularización permite mantener los parámetros del modelo controlados al agregar una penalización en la función objetivo a qué tan grande sean estos parámetros. Esto es controlado por una constante \(\lambda \in [0,1]\) que controla la cantidad de overfitting del modelo. En general, penalizar tanto los parámetros genera un riesgo de underfitting, ya que no permite que los parámetros crezcan lo suficiente, por eso se suele ajustar usando datos de validación.

Pregunta 2: Directed Graphical Models

Diseñe una DAG para un problema causal inventado por usted de al menos 5 nodos (puede basarse en el ejemplo de Eugene Charniak) usando Dagitty y considere que la DAG tenga al menos: una chain, un fork y un collider. Muestre con dagitty todas las independencias condicionales de su DAG. Explique las independencias usando las reglas de d-separación.

install.packages(dagitty)

Respuesta

Supongamos que estamos modelando la relación entre el clima, la cantidad de ejercicio, el estado de ánimo, el nivel de energía y la productividad diaria.

Clima: Representa el clima del día. Ejercicio: Indica la cantidad de ejercicio realizado. Estado de Ánimo: Representa el estado de ánimo del individuo. Energía: Nivel de energía. Productividad: Nivel de productividad diaria.

library(dagitty)

dag <- dagitty("dag {
  Clima -> Ejercicio 
  Ejercicio -> Energia
  Clima -> Estado_de_Animo 
  Productividad -> Estado_de_Animo 
  Productividad -> Energia
}")


# Visualización del DAG
plot(dag)
## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.

Chain (C -> E -> N):

El clima (C) y la energía (N) son d-separados dado el ejercicio (E). Conociendo la cantidad de ejercicio realizada, el clima no proporciona información adicional sobre la energía, y viceversa.

Fork (M -> P <- N):

El estado de ánimo (M) y el nivel de energía (N) son d-separados dado la productividad (P). Si conocemos la productividad, el estado de ánimo y la energía son condicionalmente independientes entre sí.

Collider (M -> P <- E):

El estado de ánimo (M) y el ejercicio (E) no son d-separados dado la productividad (P). En este caso, si conocemos la productividad, el estado de ánimo y el ejercicio están condicionalmente dependientes entre sí.

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.3.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(tidyr)

# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
# Análisis bayesiano
library("rethinking")
## Loading required package: cmdstanr
## This is cmdstanr version 0.6.1
## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
## - Use set_cmdstan_path() to set the path to CmdStan
## - Use install_cmdstan() to install CmdStan
## Loading required package: posterior
## Warning: package 'posterior' was built under R version 4.3.2
## This is posterior version 1.5.0
## 
## Attaching package: 'posterior'
## 
## The following objects are masked from 'package:stats':
## 
##     mad, sd, var
## 
## The following objects are masked from 'package:base':
## 
##     %in%, match
## 
## Loading required package: parallel
## rethinking (Version 2.40)
## 
## Attaching package: 'rethinking'
## 
## The following object is masked from 'package:purrr':
## 
##     map
## 
## The following object is masked from 'package:stats':
## 
##     rstudent

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
install.packages(c("mvtnorm","loo","coda","remotes","shape"), repos="https://cloud.r-project.org/",dependencies=TRUE)
remotes::install_github("rmcelreath/rethinking")

Pregunta 3: Regresión Lineal Bayesiana

El objetivo de esta pregunta es introducirlo en la aplicación de una regresión bayesiana. Con esto, buscaremos entender como calcular una regresión bayesiana en base al “motor” aproximación de Laplace, revisando como se comporta la credibilidad de sus predicciones y como la regresión lineal puede llegar a mutar a aplicar una transformación en el vector \(x\). Para responder esta pregunta centre su desarrollo solo en las clases y las funciones que nos brinda la librería rethinking.

Unos expertos en alometría deciden realizar un estudio de las alturas de unos niños en un colegio, buscando generar un regresor lineal bayesiano capaz de predecir la altura en base al peso de los alumnos. Para realizar este trabajo recopilan los datos table_height.csv, quien posee observaciones fisiológicas de 192 alumnos.

Parte I

En conocimiento los datos recopilados por los expertos, le solicitan realizar la siguiente serie de tareas:

# Cargamos los datos
data <- read.csv("table_height.csv")

hist(data$age)

hist(data$height)

hist(data$weight)

#para la creacion de los priors, encontraremos hacia donde se aproxima la altura
print('Mean')
## [1] "Mean"
print(mean(data$height)) #media para el b0
## [1] 108.3189
print('Standard deviation')
## [1] "Standard deviation"
print(sd(data$height)) #desviacion estandar para el b0
## [1] 25.74514
# Definimos el modelo de regresión bayesiana con priors
model <- quap(
  alist(
    height ~ dnorm(mu, sigma),  # Likelihood
    mu <- b0 + b1 * weight,       # Modelo lineal
    b0 ~ dnorm(108, 25),         # Prior para intercepto
    b1 ~ dnorm(0, 1),           # Prior para pendiente
    sigma ~ dunif(0, 25)        # Prior para desviación estándar
  ),
  data = data
)

# Ajustamos el modelo con Laplace approximation
precis(model, prob = 0.95)

Vemos que el promedio de b0 es 58.62cm y b1 aumenta 2.7cm por cada kilo.

# Extraemos los parámetros del modelo
posterior_samples <- extract.samples(model)

weight_seq <- seq(from=4 , to=50 , by=1 )
# MAP
map <- link(model, data = list(weight = weight_seq))

# Simulamos alturas no promedio
height_sim <- function(weight) {
  rnorm(
    n = nrow(posterior_samples),
    mean = posterior_samples$b0 + posterior_samples$b1 * weight,
    sd = posterior_samples$sigma
  )
}

simulated_heights <- sapply(weight_seq, height_sim)


# Calculamos intervalos de HPDI
height.HPDI <- apply(simulated_heights, 2, HPDI, prob = 0.95)

# Plot de los datos crudos
plot(data$weight, data$height, col = rgb(0, 0, 0, 0.5), xlab = "Weight", ylab = "Height")

# Dibujamos la línea MAP
lines(weight_seq, apply(map, 2, mean), col = "blue", lwd = 2)


# Dibujamos el intervalo HPDI para la línea
shade(apply(map, 2, HPDI, prob = 0.95), weight_seq, col = rgb(0, 0, 1, 0.2))

# Dibujamos el intervalo HPDI para alturas simuladas
shade(height.HPDI, weight_seq, col = rgb(1, 0, 0, 0.2))

# Añadimos detalles al gráfico
legend("topright", legend = c("MAP Line", "HPDI", "HPDI Simulated Heights"),
       col = c("blue", rgb(0, 0, 1, 0.2), rgb(1, 0, 0, 0.2)), lwd = c(2, 0, 0), pch = c(NA, NA, NA))

Vemos que los valores de peso vs altura siguen una tendencia de estar cerca de la linea central, los limites de la sombra encierran a la gran mayoria de los datos, y la linearidad se observa en esto Parte II

En base a los resultados obtenidos, el experto que trabaja con usted le señala que las alturas se suelen modelas con pesos logarítmicos, por lo que le sugiere añadir un logaritmo natural en el vector \(x\) que compone su modelo lineal. Realice nuevamente la regresión utilizando un intervalo del \(95\%\) sobre la media y los valores predecidos de la altura. Comente los resultados obtenidos, señalando si el modelo logra ajustar mejor los valores.

Respuesta:

Tenemos que realizar todo lo que hicimos previamente, asi que me saltare bastantes pasos:

data$weight <- log(data$weight)

hist(data$weight)

#para la creacion de los priors
print('Mean')
## [1] "Mean"
print(mean(data$weight)) #media para el b0
## [1] 2.793843
print('Standard deviation')
## [1] "Standard deviation"
print(sd(data$weight)) #desviacion estandar para el b0
## [1] 0.5010516
# Definimos el modelo de regresión bayesiana con priors
model <- quap(
  alist(
    height ~ dnorm(mu, sigma),  # Likelihood
    mu <- b0 + b1 * weight,       # Modelo lineal
    b0 ~ dnorm(108, 25),         # Prior para intercepto
    b1 ~ dnorm(20, 10),           # Prior para pendiente
    sigma ~ dunif(0, 25)        # Prior para desviación estándar
  ),
  data = data
)

posterior_samples <- extract.samples(model)

weight_seq <- seq(from=0.6 , to=log(50) , by=0.4 )
# MAP
map <- link(model, data = list(weight = weight_seq))

# Simulamos alturas no promedio
height_sim <- function(weight) {
  rnorm(
    n = nrow(posterior_samples),
    mean = posterior_samples$b0 + posterior_samples$b1 * weight,
    sd = posterior_samples$sigma
  )
}

simulated_heights <- sapply(weight_seq, height_sim)


# Calculamos intervalos de HPDI
height.HPDI <- apply(simulated_heights, 2, HPDI, prob = 0.95)

# Plot de los datos crudos
plot(data$weight, data$height, col = rgb(0, 0, 0, 0.5), xlab = "Weight", ylab = "Height")

# Dibujamos la línea MAP
lines(weight_seq, apply(map, 2, mean), col = "blue", lwd = 2)


# Dibujamos el intervalo HPDI para la línea
shade(apply(map, 2, HPDI, prob = 0.95), weight_seq, col = rgb(0, 0, 1, 0.2))

# Dibujamos el intervalo HPDI para alturas simuladas
shade(height.HPDI, weight_seq, col = rgb(1, 0, 0, 0.2))

# Añadimos detalles al gráfico
legend("topright", legend = c("MAP Line", "HPDI", "HPDI Simulated Heights"),
       col = c("blue", rgb(0, 0, 1, 0.2), rgb(1, 0, 0, 0.2)), lwd = c(2, 0, 0), pch = c(NA, NA, NA))

Podemos notar que al realizar este proceso logaritmico los resultados obtenidos son más consistentes con la regresion lineal, ya que estos se concentran más dentro de la linea y su sombra

Bonus

Esta parte es opcional y contará como un puntaje adicional si la realizan, pero si no lo hacen no habrá penalización.

Compare los resultados obtenidos con el método de Laplace con MCMC y comente los resultados. Para esto pueden utilizar ulam del paquete rethinking.

install.packages("rstan", repos = c("https://mc-stan.org/r-packages/", getOption("repos"))) #necesitarán instalar estas librerías
install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
# si no les funciona instalar cmdstan con lo anterior pueden probar con:
check_cmdstan_toolchain(fix = TRUE)
install_cmdstan()

Respuesta Aquí

Pregunta 4: MCMC

El objetivo de esta pregunta es lograr samplear, mediante la técnica de MCMC, la distribución gamma.

En general la distribución gamma se denota por \(\Gamma(\alpha,\beta)\) donde \(\alpha\) y \(\beta\) son parámetros positivos, a \(\alpha\) se le suele llamar “shape” y a \(\beta\) rate La densidad no normalizada de una distribución gamma esta dada por:

\[ f(x\mid \alpha,\beta) = \begin{cases} x^{\alpha -1}e^{-\beta x} ~ &\text{ si } x > 0\\ 0 ~&\text{si } x \leq 0 \end{cases} \]

  • Implemente el algoritmo de metropolis hasting utilizando \(q(\theta^{(t)} \mid \theta^{(t-1)}) = \mathcal{N}(\theta^{t-1},1)\), importante su función tiene que poder recibir:
    • La condición inicial \(\theta_{0}\).
    • Cantidad de repeticiones.
    • \(\alpha\) y \(\beta\).
    Escriba el algoritmo sin utilizar implementenaciones de la distribución gamma en r.

De ahora en adelante considere \(\alpha = 5\) y \(\beta = \frac{1}{5}\).

  • Considere \(\theta_{0} = 1\), grafique el histograma que obtiene para distintas cantidad de repeticiones, entre sus pruebas tiene que tener al menos una que tenga del orden de \(10^5\) repeticiones ¿Como cambia la distribución en funcion de las repeticiones?
  • Considere distintos valores de \(\theta_{0}\) y una cantidad de repeticiones grande (del orden de \(10^5\)), grafique las distribuciones que obtenga comente sus resultados ¿es importante la condición inicial del algoritmo?.
  • Utilizando \(\theta_{0}\) y cantidad de repeticiones conveniente (de acuerdo a lo que obtuvo en las partes anteriores) compare con la distribución real. Obs: En esta parte puede utilizar la distribución gamma que tiene implementado r para comparar con lo que usted realizo.

Respuesta:

Primero, definimos la función gamma_func para calcular los valores de la función de densidad de la distribución gamma, según la definición dada anteriormente \[ f(x\mid \alpha,\beta) = \begin{cases} x^{\alpha -1}e^{-\beta x} ~ &\text{ si } x > 0\\ 0 ~&\text{si } x \leq 0 \end{cases} \]

gamma_func <- function(x, shape, rate) {
  if (x > 0)
    return(x^(shape - 1) * exp(-rate * x))
  else
    return(0)
}

Luego, definimos la función metropolis la cual nos va a permitir calcular lo pedido

metropolis <- function(theta_0, n_iter, alpha, beta) {
  theta <- vector(length = n_iter + 1)
  current <- theta_0
  for (k in 1:n_iter + 1) {
    theta[k] <- current
    # Obtenemos el valor propuesto usando la distribución normal
    proposal <- rnorm(1, mean = theta[k], sd = 1)
    # Valor de la densidad gamma para este valor propuesto
    prob_prop <- gamma_func(proposal, shape = alpha, rate = beta)
    # Probabilidad para el valor actual
    prob_current <- gamma_func(current, shape = alpha, rate = beta)
    # Proporción entre ambas probabilidades
    ratio <- prob_prop / prob_current
    # Probabilidad de aceptar el nuevo valor
    prob <- min(1, ratio)
    # Decisión sobre si aceptar el nuevo valor o no
    decision <- rbinom(1, 1, prob)
    # Actualización del valor
    current <- ifelse(decision == 1, proposal, current)
  }
  # Devolvemos theta_1 a theta_(n+1), dando n valores
  return(theta[-1])
}

Ahora, graficamos el histograma para la cantidad de valores \(10^1\), \(10^2\), \(10^3\), \(10^4\), \(10^5\) y \(10^6\), considerando \(\theta_0 = 1\)

hist(metropolis(1, 1e1, 5, 0.2))

hist(metropolis(1, 1e2, 5, 0.2))

hist(metropolis(1, 1e3, 5, 0.2))

hist(metropolis(1, 1e4, 5, 0.2))

hist(metropolis(1, 1e5, 5, 0.2))

hist(metropolis(1, 1e6, 5, 0.2))

Notemos que a partir de las 1000 iteraciones, el gráfico se comienza a parecer a la distribución gamma. Mientras más iteraciones, mejor se ajusta el modelo.

Ahora cambiemos la condición inicial. Usaremos los valores \(\theta \in \{1, 5, 10, 50, 100\}\) con \(n_{iter} = 5 * 10^5\)

hist(metropolis(1, 5e5, 5, 0.2))

hist(metropolis(5, 5e5, 5, 0.2))

hist(metropolis(10, 5e5, 5, 0.2))

hist(metropolis(50, 5e5, 5, 0.2))

hist(metropolis(100, 5e5, 5, 0.2))

Notemos que en general la forma de los gráficos no cambia dependiendo de la condición inicial. El algoritmo MCMC converge igualmente a la distribución buscada. Lo que podría pasar es que, si la condición inicial es “muy mala”, el algoritmo necesite más iteraciones para converger.

Ahora, utilizando \(\theta_0 = 1\) y \(n_{iter} = 10^6\), comparemos el gráfico con la función gamma nativa de R, rgamma.

Para esto, graficaremos la función de densidad gamma usando la función rgamma. Para esto, generaremos \(10^6\) elementos aleatorios considerando \(\alpha = 5\) y \(\beta = \frac{1}{5} = 0.2\)

# Distribución no normalizada
gamma_vect <- rgamma(1e6, shape = 5, rate = 0.2)
hist(gamma_vect, xlim = range(c(0, 120)))

Luego, grafiquemos usando la función metropolis definida

# Distribución no normalizada
gamma_vect <- metropolis(1, 1e6, 5, 0.2)
hist(gamma_vect, xlim = range(c(0, 120)))

Notemos que, efectivamente, ambas gráficas son bastante similares en forma, lo que nos permite afirmar que el algoritmo metropolis se ajusta bien a la distribución gamma y por lo tanto permite estimarla correctamente.

 

A work by CC6104